rm(list=ls())

# This code replicates Figure 2 
# and accompanying Appendix tables

pacman::p_load(tidyverse, estimatr, 
               simcf, interflex,
               texreg, patchwork) 
# tidyverse


# rescale function ----
rescale01 <- function(x){
  minx <- min(x, na.rm=T)
  maxx <- max(x, na.rm=T)
  return((x-minx)/(maxx - minx))
}

# START ANALYSIS
df <- read_csv('clean_survey2.csv')

# model 1 Replication Lucid Data main findings Figure C11 ---- 
summary(lm1 <- lm_robust(policy_scale ~ sea_level_rise +
                           pid7_r + ideo5_c + age + female + college + 
                           income_60_125 + income_over_125 + income_missing + 
                           white + factor(wave), df, cluster=fips))

texreg::texreg(list(lm1),
               custom.coef.names = c('Intercept','Susceptibility','Party ID (R)', 
                                     'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White', 'Wave 2', 'Wave 3'),
               custom.model.names = c('Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               file = 'tables/appendix_table_main_model_replication_lucid.tex')

# models for identity for Figure 2 ----

# willingness to move ----
form <- policy_scale ~ sea_level_rise*move_not_at_all +
  sea_level_rise*move_not_too +
  sea_level_rise*move_somewhat +
  pid7_r + ideo5_c + age + female + college + 
  income_60_125 + income_over_125 + income_missing + 
  white + wave2 + wave3
summary(lm1c <- lm_robust(form, df, cluster=fips))
screenreg(list(lm1c))
pe <- coef(lm1c)
vc <- vcov(lm1c)
set.seed(12345)
simbetas <- MASS::mvrnorm(10000, pe, vc)

xhyp <- seq(0,1,length.out=2)
mm_low <- cfMake(form, df, nscen=2)
mm_low <- cfChange(mm_low, 'sea_level_rise', xpre=0,x=1, scen=1)
mm_low <- cfChange(mm_low, 'move_not_at_all', xpre=1,x=1, scen=1)
mm_low <- cfChange(mm_low, 'move_not_too', xpre=0, x=0, scen=1)
mm_low <- cfChange(mm_low, 'move_somewhat', xpre=0, x=0, scen=1)

mm_low <- cfChange(mm_low, 'sea_level_rise', xpre=0,x=1, scen=2)
mm_low <- cfChange(mm_low, 'move_not_at_all', xpre=0, x=0, scen=2)
mm_low <- cfChange(mm_low, 'move_not_too', xpre=0, x=0, scen=2)
mm_low <- cfChange(mm_low, 'move_somewhat', xpre=0, x=0, scen=2)

willing_move <- linearsimfd(mm_low, simbetas,ci=0.9) %>% 
  as_tibble() %>% 
  mutate(moderator = c('Least Willing','Most Willing'))

# ties to community ----
form <- policy_scale ~ sea_level_rise*comm_ties_strong +
  sea_level_rise*comm_ties_sw_strong + 
  sea_level_rise*comm_ties_nv_strong +
  pid7_r + ideo5_c + age + female + college + 
  income_60_125 + income_over_125 + income_missing + 
  white + wave2 + wave3
summary(lm3c <- lm_robust(form, df, cluster=fips))

pe <- coef(lm3c)
vc <- vcov(lm3c)
set.seed(12345)
simbetas <- MASS::mvrnorm(10000, pe, vc)

xhyp <- seq(0,1,length.out=2)
mm_low <- cfMake(form, df, nscen=2)
mm_low <- cfChange(mm_low, 'sea_level_rise', xpre=0,x=1, scen=1)
mm_low <- cfChange(mm_low, 'comm_ties_strong', xpre=1,x=1, scen=1)
mm_low <- cfChange(mm_low, 'comm_ties_sw_strong', xpre=0, x=0, scen=1)
mm_low <- cfChange(mm_low, 'comm_ties_nv_strong', xpre=0, x=0, scen=1)

mm_low <- cfChange(mm_low, 'sea_level_rise', xpre=0,x=1, scen=2)
mm_low <- cfChange(mm_low, 'comm_ties_strong', xpre=0, x=0, scen=2)
mm_low <- cfChange(mm_low, 'comm_ties_sw_strong', xpre=0, x=0, scen=2)
mm_low <- cfChange(mm_low, 'comm_ties_nv_strong', xpre=0, x=0, scen=2)

comm_ties <- linearsimfd(mm_low, simbetas,ci=0.9) %>% 
  as_tibble() %>% 
  mutate(moderator = c('Strongest','Weakest'))

# place based identity ----
out.if <- interflex(Y = "policy_scale", 
                    D = "sea_level_rise", 
                    X = "place_identity_2",
                    Z=c('pid7_r','ideo5_c','age','female','college',
                        'income_60_125','income_over_125','income_missing',
                        'white', 'wave2','wave3'),
                    cl='fips',
                    data = as.data.frame(df), 
                    estimator = "binning", 
                    vcov.type = "cluster", 
                    main = "Marginal Effects",
                    nbins = 3,
                    na.rm=T)


# extract model
summary(lm4c <- out.if$model.binning)
lm4cSE <- sqrt(diag(vcovCluster(lm4c, out.if$model.binning$data$fips)))
lm4cPV <- lmtest::coeftest(lm4c, vcovCL(lm4c, cluster = out.if$model.binning$data$fips))[,4]

# save interflex figure
interflex_identity <- out.if$est.bin$`sea_level_rise=0 (50%)` %>%
  as_tibble() %>%
  ggplot(aes(x=x0, y=coef, ymin=CI.lower,ymax=CI.upper)) +
  geom_ribbon(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, ymin=`lower CI(95%)`,ymax=`upper CI(95%)`), 
              inherit.aes = F, fill='grey80') +
  geom_line(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, y=ME), 
            inherit.aes = F, color='grey50') +
  geom_hline(yintercept=0, color='red',linetype=2) +
  geom_pointrange() +
  theme_bw() +
  labs(x='Moderator (Place-Based Identity)',
       y='Marginal Effect SLR on\nPolicy Attitudes')

# save binning marginal effects
place_id <- data.frame(
  pe = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,2],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,2]), 
  se = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,3],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,3]),
  moderator = c('Lowest','Highest')
)

# 90% CIs
place_id$lower <- place_id$pe - qnorm(0.95) * place_id$se
place_id$upper <- place_id$pe + qnorm(0.95) * place_id$se

# plot for identity items ----
identity_plot <- bind_rows(
  willing_move %>% mutate(mod_name = 'Willingness to Move'),
  comm_ties %>% mutate(mod_name = 'Community Ties'),
  place_id %>% mutate(mod_name = 'Place-Based Identity')
) %>%
  ggplot(aes(x=moderator, y=pe, ymin=lower, ymax=upper)) + 
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_errorbar(width=0.2) +
  geom_point(shape=21, size=3, fill='white') +
  facet_wrap(~mod_name, scales='free', ncol=1) +
  coord_flip() +
  scale_y_continuous(limits=c(-0.12,0.25)) +
  theme_bw() +
  labs(x='', y='Change in Support Climate Mitigation')
identity_plot # will save later

# tables C22 regressions for identity moderators -----
texreg::texreg(list(lm1c, lm3c),
               custom.model.names = c('Support Policy','Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               reorder.coef=c(1,2,6:16,3:5,17:25),
               custom.coef.names = c('Intercept','Susceptibility',
                                     'Unwill Move 4','Unwill Move 3','Unwill Move 2',
                                     'Party ID (R)', 'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White', 'Wave 2', 'Wave 3',
                                     "Susceptibility * Unwill Move 4", "Susceptibility * Unwill Move 3", "Susceptibility * Unwill Move 2",
                                     "Community Ties 4", "Community Ties 3", "Community Ties 2",
                                     "Susceptibility * Community Ties 4", "Susceptibility * Community Ties 3", "Susceptibility * Community Ties 2"),
               file = 'tables/appendix_table_identities_1.tex')

texreg::texreg(list(lm4c),
                  custom.model.names = c('Support Policy'),
                  stars = c(0.1,0.05,0.01), include.ci=F,
                  override.se = list(lm4cSE),
                  override.pvalues = list(lm4cPV),
                  custom.coef.names = c('G1','GX1','G2','GX2','G3','GX3',
                                        'DG1','DGX1','DG2','DGX2','DG3','DGX3',
                      'Party ID (R)', 'Conservative', 'Age','Female','College',
                      'Income 60-125k','Income Over 125k','Income Missing',
                      'White', 'Wave 2', 'Wave 3'),
                  file = 'tables/appendix_table_identities_2.tex')


# economic self-interest -----

# homeowner ----
form <- policy_scale ~ sea_level_rise*homeowner +
  pid7_r + ideo5_c + age + female + college + 
  income_60_125 + income_over_125 + income_missing + 
  white + wave2 + wave3
summary(lm1 <- lm_robust(form, df, cluster=fips))

pe <- coef(lm1)
vc <- vcov(lm1)
set.seed(12345)
simbetas <- MASS::mvrnorm(10000, pe, vc)

xhyp <- seq(0,1,length.out=2)
mm_low <- cfMake(form, df, nscen=2)
mm_low <- cfChange(mm_low, 'sea_level_rise', xpre=0,x=1, scen=1)
mm_low <- cfChange(mm_low, 'homeowner', xpre=1,x=1, scen=1)
mm_low <- cfChange(mm_low, 'sea_level_rise', xpre=0,x=1, scen=2)
mm_low <- cfChange(mm_low, 'homeowner', xpre=0, x=0, scen=2)

homeownership <- linearsimfd(mm_low, simbetas,ci=0.9) %>% 
  as_tibble() %>% 
  mutate(moderator = c('Homeowner','Non-Homeowner'))

# insurance models -----
out.if <- interflex(Y = "policy_scale", 
                    D = "sea_level_rise", 
                    X = "cost_insurance_log",
                    Z=c('pid7_r','ideo5_c','age','female','college',
                        'income_60_125','income_over_125','income_missing',
                        'white', 'wave2'),
                    cl='fips',
                    data = as.data.frame(df[df$homeowner==1,]), 
                    estimator = "binning", 
                    vcov.type = "cluster", 
                    main = "Marginal Effects",
                    nbins = 3,
                    treat.type = 'continuous',
                    na.rm=T)

# save interflex figure
interflex_homeinsurance <- out.if$est.bin$`sea_level_rise=0 (50%)` %>%
  as_tibble() %>%
  ggplot(aes(x=x0, y=coef, ymin=CI.lower,ymax=CI.upper)) +
  geom_ribbon(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, ymin=`lower CI(95%)`,ymax=`upper CI(95%)`), 
              inherit.aes = F, fill='grey80') +
  geom_line(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, y=ME), 
            inherit.aes = F, color='grey50') +
  geom_hline(yintercept=0, color='red',linetype=2) +
  geom_pointrange() +
  theme_bw() +
  labs(x='Moderator (Home Insurance)',
       y='Marginal Effect SLR on\nPolicy Attitudes')

# extract and save regression
summary(lm2 <- out.if$model.binning)
lm2SE <- sqrt(diag(vcovCluster(lm2, out.if$model.binning$data$fips)))
lm2PV <- lmtest::coeftest(lm2, vcovCL(lm2, cluster = out.if$model.binning$data$fips))[,4]

# save estimates for figure
insurance <- data.frame(
  pe = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,2],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,2]), 
  se = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,3],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,3]),
  moderator = c('Lowest','Highest')
)
# 90% CI
insurance$lower <- insurance$pe - qnorm(0.95) * insurance$se
insurance$upper <- insurance$pe + qnorm(0.95) * insurance$se



# income ----
out.if <- interflex(Y = "policy_scale", 
                     D = "sea_level_rise", 
                     X = "hhi_cont",
                     Z=c('pid7_r','ideo5_c','age','female','college',
                         'white','wave2','wave3'),
                     data = as.data.frame(df), 
                    cl='fips',
                     estimator = "binning", 
                     vcov.type = "cluster", 
                     main = "Marginal Effects",
                     na.rm=T)

# extract binning model
summary(lm3 <- out.if$model.binning)
lm3SE <- sqrt(diag(vcovCluster(lm3, out.if$model.binning$data$fips)))
lm3PV <- lmtest::coeftest(lm3, vcovCL(lm3, cluster = out.if$model.binning$data$fips))[,4]

# save estimates
income <- data.frame(
  pe = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,2],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,2]), 
  se = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,3],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,3]),
  moderator = c('Lowest','Highest')
)

# 90% CI
income$lower <- income$pe - qnorm(0.95) * income$se
income$upper <- income$pe + qnorm(0.95) * income$se

# tourism -----
df$n_hotels_log <- rescale01(log(df$n_hotels+1))
out.if <- interflex(Y = "policy_scale", 
                    D = "sea_level_rise", 
                    X = "n_hotels_log",
                    Z=c('pid7_r','ideo5_c','age','female','college',
                        'income_60_125','income_over_125','income_missing',
                        'white', 'wave2','wave3'),
                    data = as.data.frame(df), 
                    estimator = "binning", 
                    cl='fips',
                    vcov.type = "cluster", 
                    main = "Marginal Effects",
                    nbins = 3,
                    na.rm=T)

# save interflex figure
interflex_hotels <- out.if$est.bin$`sea_level_rise=0 (50%)` %>%
  as_tibble() %>%
  ggplot(aes(x=x0, y=coef, ymin=CI.lower,ymax=CI.upper)) +
  geom_ribbon(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, ymin=`lower CI(95%)`,ymax=`upper CI(95%)`), 
              inherit.aes = F, fill='grey80') +
  geom_line(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, y=ME), 
            inherit.aes = F, color='grey50') +
  geom_hline(yintercept=0, color='red',linetype=2) +
  geom_pointrange() +
  theme_bw() +
  labs(x='Moderator (Number Hotels)',
       y='Marginal Effect SLR on\nPolicy Attitudes')

# extract regression table
summary(lm4 <- out.if$model.binning)
lm4SE <- sqrt(diag(vcovCluster(lm4, out.if$model.binning$data$fips)))
lm4PV <- lmtest::coeftest(lm4, vcovCL(lm4, cluster = out.if$model.binning$data$fips))[,4]

# save estimates
hotels <- data.frame(
  pe = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,2],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,2]), 
  se = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,3],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,3]),
  moderator = c('Fewest','Most')
)
# 90% CI
hotels$lower <- hotels$pe - qnorm(0.95) * hotels$se
hotels$upper <- hotels$pe + qnorm(0.95) * hotels$se


# payrolls -----
hist(log(df$payroll_1000+1))
df$payroll_1000_log <- rescale01(log(df$payroll_1000+1))
out.if <- interflex(Y = "policy_scale", 
                    D = "sea_level_rise", 
                    X = "payroll_1000_log",
                    Z=c('pid7_r','ideo5_c','age','female','college',
                        'income_60_125','income_over_125','income_missing',
                        'white', 'wave2','wave3'),
                    data = as.data.frame(df), 
                    estimator = "binning", 
                    cl='fips',
                    vcov.type = "cluster", 
                    main = "Marginal Effects",
                    nbins = 3,
                    na.rm=T)

# save interflex figure
interflex_payrolls <- out.if$est.bin$`sea_level_rise=0 (50%)` %>%
  as_tibble() %>%
  ggplot(aes(x=x0, y=coef, ymin=CI.lower,ymax=CI.upper)) +
  geom_ribbon(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, ymin=`lower CI(95%)`,ymax=`upper CI(95%)`), 
              inherit.aes = F, fill='grey80') +
  geom_line(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, y=ME), 
            inherit.aes = F, color='grey50') +
  geom_hline(yintercept=0, color='red',linetype=2) +
  geom_pointrange() +
  theme_bw() +
  labs(x='Moderator (Hotel Payrolls)',
       y='Marginal Effect SLR on\nPolicy Attitudes')

# extract regression table
summary(lm5 <- out.if$model.binning)
lm5SE <- sqrt(diag(vcovCluster(lm5, out.if$model.binning$data$fips)))
lm5PV <- lmtest::coeftest(lm5, vcovCL(lm5, cluster = out.if$model.binning$data$fips))[,4]

# save estimates
payrolls <- data.frame(
  pe = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,2],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,2]), 
  se = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,3],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,3]),
  moderator = c('Lowest','Highest')
)

# 90% CI
payrolls$lower <- payrolls$pe - qnorm(0.95) * payrolls$se
payrolls$upper <- payrolls$pe + qnorm(0.95) * payrolls$se

# Home Values  -----
df$zvhi_allhomes_2020_log <- log(df$zvhi_allhomes_2020+1)
out.if <- interflex(Y = "policy_scale", 
                    D = "sea_level_rise", 
                    X = "zvhi_allhomes_2020_log",
                    Z=c('pid7_r','ideo5_c','age','female','college',
                        'income_60_125','income_over_125','income_missing',
                        'white','wave2','wave3'),
                    data = as.data.frame(df), 
                    cl='fips',
                    estimator = "binning", 
                    vcov.type = "cluster", 
                    main = "Marginal Effects",nbins = 3,
                    na.rm=T)

# save interflex object
interflex_homevalues <- out.if$est.bin$`sea_level_rise=0 (50%)` %>%
  as_tibble() %>%
  ggplot(aes(x=x0, y=coef, ymin=CI.lower,ymax=CI.upper)) +
  geom_ribbon(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, ymin=`lower CI(95%)`,ymax=`upper CI(95%)`), 
              inherit.aes = F, fill='grey80') +
  geom_line(data=out.if$est.lin$`sea_level_rise=0 (50%)`, aes(x=X, y=ME), 
            inherit.aes = F, color='grey50') +
  geom_hline(yintercept=0, color='red',linetype=2) +
  geom_pointrange() +
  theme_bw() +
  labs(x='Moderator (Home Values)',
       y='Marginal Effect SLR on\nPolicy Attitudes')

# save regression table
summary(lm6 <- out.if$model.binning)
lm6SE <- sqrt(diag(vcovCluster(lm6, out.if$model.binning$data$fips)))
lm6PV <- lmtest::coeftest(lm6, vcovCL(lm6, cluster = out.if$model.binning$data$fips))[,4]

# extract estimates
homevals <- data.frame(
  pe = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,2],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,2]), 
  se = c(out.if$est.bin$`sea_level_rise=0 (50%)`[1,3],
         out.if$est.bin$`sea_level_rise=0 (50%)`[3,3]),
  moderator = c('Lowest','Highest')
)

# 90% CI
homevals$lower <- homevals$pe - qnorm(0.95) * homevals$se
homevals$upper <- homevals$pe + qnorm(0.95) * homevals$se

# make econ plot -----
econ_plot <- bind_rows(
  homeownership %>% mutate(mod_name = 'Homeownership'),
  homevals %>% mutate(mod_name = 'Home Values'),
  insurance %>% mutate(mod_name = 'Home Insurance Cost'),
  income %>% mutate(mod_name = 'Income'),
  hotels %>% mutate(mod_name = 'Number Hotels'),
  payrolls %>% mutate(mod_name = 'Hotel Payrolls')
) %>%
  mutate(mod_name = factor(mod_name, levels=c(
    'Homeownership','Home Values','Home Insurance Cost','Income','Number Hotels','Hotel Payrolls'
  ))) %>%
  mutate(moderator = factor(moderator, levels=c(
    'Homeowner','Non-Homeowner',
    'Highest Cost','Lowest Cost',
    'Highest','Lowest',
    'Most','Fewest'
  ))) %>%
  ggplot(aes(x=moderator, y=pe, ymin=lower, ymax=upper)) + 
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_errorbar(width=0.2) +
  geom_point(shape=21, size=3, fill='white') +
  facet_wrap(~mod_name, scales='free', ncol=1) +
  coord_flip() +
  scale_y_continuous(limits=c(-0.1,0.25)) +
  theme_bw() +
  labs(x='', y='Change in Support Climate Mitigation')

# save econ tables ----
texreg::texreg(list(lm1),
               custom.model.names = c('Support Policy'),
               stars = c(0.1,0.05,0.01), include.ci=F,
               reorder.coef=c(1,2,3,15,4:14),
               custom.coef.names = c('Intercept','Susceptibility',
                                     'Homeowner',
                                     'Party ID (R)', 'Conservative', 'Age','Female','College',
                                     'Income 60-125k','Income Over 125k','Income Missing',
                                     'White', 'Wave 2', 'Wave 3',
                                     "Susceptibility * Homeowner"),
               file = 'tables/appendix_table_economics_1.tex')

texreg::texreg(list(lm2,lm3,lm4,lm5,lm6),
                  custom.model.names = c('Policy (Insurance)','Policy (Income)','Policy (Hotels)',
                  'Policy (Payrolls)','Policy (Home Values)'),
                  override.se = list(lm2SE,lm3SE,lm4SE,lm5SE,lm6SE),
                  override.pvalues = list(lm2PV,lm3PV,lm4PV,lm5PV,lm6PV),
                  stars = c(0.1,0.05,0.01), include.ci=F,
                  #reorder.coef=c(1,2,4:14,3,15:35),
                  custom.coef.names = c('G1','GX1','G2','GX2','G3','GX3',
                                        'DG1','DGX1','DG2','DGX2','DG3','DGX3',
                                        'Party ID (R)', 'Conservative', 'Age','Female','College',
                                        'Income 60-125k','Income Over 125k','Income Missing',
                                        'White', 'Wave 2', 'Wave 3'),
                  file = 'tables/appendix_table_economics2.tex')

#ggsave(econ_plot, width=4, height=5, file='figures/figure_econ_mod.png')
#ggsave(identity_plot, width=4, height=5, file='figures/figure_ident_mod.png')

econ_plot + identity_plot + plot_layout(axis_titles = "collect")
ggsave(width=9, height=5.5, file='figures/figure_moderators.pdf')


# write out interflex for R1 ----

interflex_identity + 
  interflex_homeinsurance +
  interflex_homevalues +
  interflex_hotels +
  interflex_payrolls 
ggsave(width=9, height=5, 
       filename = 'figures/interflex_estimates_for_appendix.pdf')

# check stat sign coefficients for figure 2 -----

t_stats_diff <- function(x){
  diff <- x[,'pe'][1] - x[,'pe'][2]
  sumSEsq <- (x[,'se'][1])^2 +  (x[,'se'][2])^2
  t.stat <- diff / sqrt(sumSEsq)
  p_value = 2*pt(abs(t.stat), df = Inf,lower.tail = FALSE)
  return(data.frame(
    diff = diff, 
    p = p_value))
} 

# econ scales 
homeownership_new <- homeownership %>% mutate(se = (pe - lower)/1.64) %>% dplyr::select(pe, se)
t_stats_diff(as.data.frame(homeownership_new))

t_stats_diff(homevals)
t_stats_diff(insurance)
t_stats_diff(income)
t_stats_diff(hotels)
t_stats_diff(payrolls)

# identity scales
willing_move_new <- willing_move %>% mutate(se = (pe - lower)/1.64) %>% dplyr::select(pe, se)
comm_ties <- comm_ties %>% mutate(se = (pe - lower)/1.64)

t_stats_diff(as.data.frame(willing_move_new))
t_stats_diff(as.data.frame(comm_ties))
t_stats_diff(as.data.frame(place_id))

